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Lattice dynamical methods used to predict phase transformations in crystals typically deal with 
harmonic phonon spectra and are therefore not applicable in important situations where one of the 
competing crystal structures is unstable in the harmonic approximation, such as the bcc structure 
involved in the hep to bcc martensitic phase transformation in Ti, Zr and Hf. Here we present an 
expression for the free energy that does not suffer from such shortcomings, and we show by self 
consistent ab initio lattice dynamical calculations (SCAILD), that the critical temperature for the 
hep to bcc phase transformation in Ti, Zr and Hf, can be effectively calculated from the free energy 
difference between the two phases. This opens up the possibility to study quantitatively, from first 
principles theory, temperature induced phase transitions. 



Martensitic phase transformations are common, both 
in alloys frequently used in industry, such as shape mem- 
ory alloys p] , and in the the elemental group 3 to 4 tran- 
sition metals [2] , not to mention martensitic transforma- 
tion in iron and iron-based alloys, a crucial phenomenon 
for metallurgy [3 . Thus there exists a substantial inter- 
est both from an industrial, applied and an academic 
point of view to develop accurate and effective meth- 
ods to understand and even predict martensitic phase- 
transformations. 

The hep to bcc (or a to /3) transition in Ti, Zr and 
Hf is a martensitic phase transformation that has been 
thoroughly investigated both from an experimental [21 H] 
and theoretical [SHE] perspective. Recently Hennings et 
al developed and used a classical potential of the modi- 
fied embedded atom method (MEAM) [5] to accurately 
reproduce the phase boundary between the hep and bcc 
structure in Ti. However, there is up to this date a lack 
of first principles theoretical studies made of the marten- 
sitic hep to bcc phase-transformation in Ti, Zr and Hf. 
The problem is that anharmonic effects in lattice dynam- 
ics |10| are of crucial importance for finite-temperature 
structural phase transitions, and their quantitative first- 
principle treatment is a real challenge. 

A straightforward calculation using DFT molecular dy- 
namics (DFT-MD) 11 should in principle be able to re- 
produce the bcc to hep phase transformation in Ti and 
similar materials, since DFT-MD implicitly include an- 
harmonic effects. However, DFT-MD is a computation- 
ally very demanding task which makes its use problem- 
atic. Instead we will here exploit the method of self con- 
sistent ab initio lattice dynamical calculations (SCAILD) 
|12j . Here, we further develop this method in order 
to be able to calculate thermodynamic properties, such 
as structural free energy difference (before we were re- 
stricted by the calculations of temperature-dependent 



phonon frequencies only [T2HIB])- Since the SCAILD 
scheme is a constrained sampling method, in that it only 
samples the lattice dynamical phase-space along the nor- 
mal mode directions of commensurate phonons [121 113) , 
the SCAILD calculations are much faster and, thus, much 
more practical than the corresponding DFT-MD calcu- 
lations. 

Thus, we propose here a new expression of the free en- 
ergy of the SCAILD scheme, and we show that from the 
atomic configurations and the phonon density of states 
produced by the SCAILD calculations, an accurate mea- 
sure of the free energy for the different phases can be 
obtained. 

In order to properly describe temperature driven phase 
transformations in general, one must include the interac- 
tion between phonons |10) . As a result, phonon frequen- 
cies turn out to be temperature dependent which we ex- 
plore numerically in this study by means of the SCAILD 
method [IMS]- 

The SCAILD method is based on the calculation of 
Hellman-Feynman forces of atoms in a supercell. The 
method can be viewed as an extension of the frozen 
phonon method |17j . in which all phonons with wave 
vectors q commensurate with the supercell are excited 
together in the same cell by displacing atoms situated at 
the undistorted positions R + b^, according to R + b CT —> 
R + h a + Ur„ , where the displacements are given by 

U ^ = -7=E^<s eiq(R+bCT) - C 1 ) 

Here R represent the N Bravais lattice sites of the super- 
cell, b CT the position of atom a relative to this site, e^ s 
are the phonon eigenvectors corresponding to the phonon 
mode, s, and the mode amplitude A° is calculated from 
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the different phonon frequencies w qs through 




(2) 



where n qs = n(-^p), with n(x) = l/(e x — 1), are the 
phonon occupational numbers, M a the atomic masses 
and T is the temperature of the system. The phonon fre- 
quencies, w qs , are defined through the variational deriva- 
tive of the total energy with respect to the occupation 
numbers 
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obtained through the Fourier transform Fq of the forces 
acting on the atoms in the supercell. 

Due to the simultaneous presence of all the commen- 
surate phonons in the same force calculation, the inter- 
action between different lattice vibrations are taken into 
account and the phonon frequencies given by Eq. ([3| are 
thus renormalized by the very same interaction. 

By alternating between calculating the forces on the 
displaced atoms and calculating new phonon frequencies 
and new displacements through Eqs. ([l])-([3| the phonon 
frequencies are calculated in a self consistent manner. For 
more details on the SCAILD method we refer to Refs. 
IT2HI41 It should be mentioned that we do not consider 
here the phonon decay processes (see, e.g., Ref. [TBI and 
references therein). 

The free energy as a function of volume, V, and tem- 
perature for the bcc and hep structures can be calculated 
through the expression 

F(T,V) = U (V)+F ph (V,T) + F el (V,T), (4) 

where Uq is the static ground state energy of the respec- 
tive structures at T = K (i.e without any phonons 
excited and temperature excitations of electronic states 
), F p h is the free energy of the phonons and F e i is the 
free energy of the electrons. The temperature dependent 
parts of the free energy can be found as 

F ph (V,T) + F el (V,T) = 

I- ]T AF*({U R },V,T) + h B T-TS ph (V,T). (5) 
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Here AF* is the change in free energy relative to the 
ground state energy Uq, caused by the phonon induced 
atomic displacements described by Eq. ([I]), and thermal 
excitations of the electronic states. The sum on the right- 
hand part of Eq. ^ is over the different atomic config- 
urations, {Ur}, generated throughout the SCAILD self 
consistent run. Since A_F* are calculated at atomic con- 
figurations accommodating the different frozen phonon 



superposition of Eq. , AF* not only contains the finite 
temperature contribution to the electronic free energy for 
a given atomic configuration, but also, the potential en- 
ergy provided by the frozen lattice waves, i.e the potential 
energy of the phonons at a particular phonon superposi- 
tion [15] . The phonon kinetic energy is given by 3fcsT/2 
per atom which means that atomic motion is considered 
as classical; typically, temperatures of structural phase 
transformations are higher than the Debye temperature, 
thus, this approximation is well justified. 

In practice, the sum of the finite temperature electron 
free energy and phonon potential energy, AF*, was ob- 
tained by calculating the total free energy of the corre- 
sponding atomic configuration using a Fermi-Dirac tem- 
perature smearing of the Kohn-Sham occupational num- 
bers [5D] , and then subtracting the static potential energy 
Uq (Here Uq is calculated with the tetrahedron method to 
provide a good reference to the temperature excited elec- 
tronic states). The number of configurations, TVj, used 
for each volume and temperature was typically 400. 

Another problem is how to calculate the phonon en- 
tropy. We will assume that it depends on the phonon 
occupation numbers n qs in the same way as for nonin- 
teracting bosons: 
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This is the only entropy expression consistent with Eq. 
@ and {3|, which can be proved in the exact same man- 
ner as was done for fermions in the Landau theory of a 
normal Fermi liquid [21]. Thus the SCAILD scheme, to- 
gether with a free energy defined through Eqs. ^ and 
([6]), constitutes nothing but a theory of a "normal Bose 
liquid" . Expression Q can be written in terms of the 
phonon density of states, <?(w), produced by a converged 
SCAILD calculation, and is given by 

TS ph (V,T) = 

du)g(u,V,T)nu\n(-r^=) - ln fl - 
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(l-e *bt^ .(7) 



Here the phonon frequencies used to calculate the phonon 
density of states, g(u)), are the normal mode configura- 
tional mean values 
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It should be stressed that the partitioning of the free 
energy through Eqs. Q, ^ and |7| has been chosen 
to maximize both the accuracy of the phonon potential 
energy, which in the form of Eqn.([5| take into account 
anharmonicity up to infinite order, and the phonon en- 
tropy, which in the form given by Eqn.([7| is accurate to 
leading order in anharmonic perturbation theory [22] . 
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The phonon density of states and the corresponding 
free energies for the hep and bec structures were cal- 
culated at up to five different temperatures, and at each 
temperature, SCAILD calculations were performed at up 
to five different volumes. As an example of the typical 
data obtained from the calculations, the resulting Tita- 
nium free energies obtained through Eq. are shown 
in Fig. g 

As regards the other computational details of the 
force calculation we used the VASP package [2U], within 
the generalized gradient approximation (GGA). The 
projector-augmented wave (PAW) potentials required en- 
ergy cutoffs of 232 eV. The Ti(4s,3d), Zr(4s,4p,5s,4d) and 
Hf(6s,5d) levels were treated as valence electrons. The k- 
point mesh was a 5 x 5 x 5 Monkhorst-Pack grid in the bec 
phase calculations. In the hep phase calculations 6x6x6 
gamma centered mesh was used. In order to include the 
electron entropy in the calculations, Fermi-Dirac temper- 
ature smearing were applied to the Kohn-Sham occupa- 
tional numbers. The bec and hep supercells used were 
obtained by increasing the bec primitive cell 4 times and 
the hep primitive cell 3 times, along the respective bec 
and hep primitive lattice vectors. The sizes of the super- 
cells where chosen such that to ensure a sufficient decay 
of the interatomic force constant within the supercell, 
permitting a proper sampling of the lattice dynamical 
phase space [23] . For the calculations of the static poten- 
tial energy, Uq in Eqn.Q, the all-electron full-potential 
linearized augmented-plane wave (FP-LAPW) package 
ELK [23] was used within the GGA approximation. This 
was found necessary to ensure a high accuracy of the zero 
temperature part of the bec - hep energy difference. An 
energy cutoff of 270 eV together with a 24 x 24 x 24 k- 
point mesh and a 24 x 24 x 15 k-point mesh were used for 
the bec and hep structures, respectively, in these calcula- 
tions the Methfessel-Paxton integration scheme was used 
with a 0.2 eV smearing of the Kohn-Sham eigenvalues. 

For each temperature the free energy obtained through 
Eq. ([5]) was fitted to a first order polynomial in V (for 
typical data see Fig. [TJ. Then by using these first order 
fits together with Eq. Q, the total free energy at each 
temperature was obtained through minimization with re- 
spect to volume. In Fig. [5Ja), (c) and [3] (a), the min- 
imized free energy at each temperature is displayed for 
the bec and hep structure in Ti, Zr and Hf, respectively. 
In Fig. [2jL), (d) and [3] (b) the free energy difference 
between the structures is displayed for Ti, Zr and Hf, 
respectively. 

The temperature driven hep to bec phase- 
transformation in Ti can bee seen to occur at T ~ 1100 
K, which is reasonably close to the experimentally 
observed phase transition temperature of 1155 K [26] . 
whereas in Zr and Hf, the transition is predicted to 
occur at T~ 920 K and T-1660 K, respectively. The 
theoretical estimates of the transition temperature in Zr 
and Hf, with their respective experimental data of 1135 
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FIG. 1: (Color online) The calculated free energy at different 
volumes and temperatures for Ti in the bec phase (a), and 
the hep phase (b). 
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FIG. 2: (Color online) In (a) and (c) the calculated free energy 
of the bec and hep phase, in Ti and Zr, respectively. In (b) and 
(d) the calculated free energy difference, AF = Fhcp — Ftcc, 
in Ti and Hf, respectively 



K and 2015 [2J, does not correspond to experimental 
data as well as in the case of Ti, however, given that first 
principles calculations (at T = 0), with current exchange 
and correlations functionals, have a problem in resolving 
energy differences between different crystallographic 
phases better than ~10 meV/atom, one may not expect 
from any first principles based theory, like the SCAILD 
method used here, to reproduce temperature induced 
phase transitions with an accuracy better than a few 
hundred Kelvin. 

It should be noted here that the presently used free en- 
ergy expression differs significantly from the previously 
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FIG. 3: (Color online) In (a) the calculated free energy in Hf 
for the bcc and hep phase represented by red empty squares 
and filled black diamonds, respectively. In (b) the calculated 
free energy difference in Hf, AF = F hcp — F bcc . 



used expression from quasi-harmonic theory (e.g. used 
in Ref. 13J ) in that it takes into account the anharmonic 
contributions to the lattice dynamical potential energy 
without projecting it down on a quasi-harmonic descrip- 
tion as in Ref. [13] ■ In fact, using the old free energy 
expression together with the first principles inter atomic 
forces calculated in this work, results in the hep phase 
in Ti being > 90 meV lower in free energy compared to 
the free energy of the bcc phase throughout the entire 
temperature intervall 800 K < T < 1200JC. 

In summary, we introduce here an expression for the 
free energy which is applicable also for highly anharmonic 
crystals. The free energy expression can readily be in- 
terfaced with lattice dynamic methods, for instance the 
SCAILD technique, and we show that the temperature 
induced hep — > bcc transition of Ti, Zr and Hf can be 
reproduced by theory. Theory puts the transition tem- 
perature in Ti, Zr and Hf to 1100 K, 920 K and 1660 K, 
respectively. 
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